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The possible emergence of a spin liquid phase in the half-filled Hubbard model on the honeycomb 
lattice; a simple model of graphene, is studied using the variational cluster approximation. We 
found that the critical interaction strength of a magnetic transition is slightly lower than that of 
the non-magnetic metal-to-insulator transition and the magnetic order parameter is already non- 
negligible at the latter transition point. Thus a semi-metallic state becomes a magnetic insulator 
as the interaction strength increases and a spin liquid state, characterized by an insulating state 
without long range order, or a state very close to that is not realized in this system. 
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Introduction .— A spin liquid state, which is insulating 
without showing any long range order, has attracted a 
lot of interests. An insulator without long range order is 
fairly common among band insulators, where the effect 
of the electron correlations are very small. An important 
difference of a spin liquid state from band insulators is 
that a spin liquid state is insulating due to the effect of 
electron correlations. In many cases, the electron corre¬ 
lations derive a system to a state with long range order, 
e.g., a magnetic state, therefore it is difficult to realize a 
spin liquid state in real materials or realistic theoretical 
models. One possibility to realize a spin liquid state is 
to consider systems with strong frustrations. In fact this 
possibility is realized in the charge organic transfer salts 
k-(BEDT-TTF) 2X [T]. Another possibility is to consider 
systems which do not contain many electrons near the 
Fermi level. In this situation, the system might become 
insulating as the degree of the electron correlations in¬ 
creases, before it becomes an ordered state due to the 
cooperative correlations of the electrons near the Fermi 
level. 

As an example of the latter possibility, the half-filled 
Hubbard model on the honeycomb lattice (See Fig. 
has been extensively studied[2lH2). The Hamiltonian of 
this model is given by 

H (I) 

<b> * 

where annihilates (creates) an electron with spin a 
on the site i, Uicr = ; U is the on-site Coulomb 

repulsion, pL is the chemical potential, and {ij) denotes 
the nearest-neighbor pairs on the lattice. This is a simple 
model describing a graphene and there are also other ma¬ 
terials exhibiting the hexagonal structures. When [7 = 0 
the band structure exhibits Dirac cones, and at half¬ 
filling the Fermi level crosses the conical vertex of the 
cones. So the density of state is zero only at the Fermi 
level and a semi-metal is realized. (See. Fig.[^(c).) Thus 
this model hosts the second possibility. 

So far the results of the preceding studies [JHE] seem 
to be rather controversial. Therefore we study this 



FIG. 1. (Color online) The honeycomb lattice. In the 
magnetic state the spins orient in the opposite directions on 
the filled and unfilled circles. The 10-site cluster circled by 
the dash-dotted line and 16-site cluster circled by the dotted 
line are used in our analysis. 

model again using the variational cluster approximation 
(VGA). VGA is classified as a kind of quantum cluster 
methods [T3], and among the previous studies the 

quantum cluster methods, including VGA, are used in 
Refs. [6HI2]. Therefore to clarify the significance of our 
study, we first briefly explain the quantum cluster meth¬ 
ods and discuss the controversial situations up to now, 
mentioning the results of previous studies. 

In the quantum cluster methods, we divide the original 
infinite lattice into identical finite clusters and approxi¬ 
mately compute physical quantities of the infinite system 
by solving a Hamiltonian defined on this finite size clus¬ 
ter. This Hamiltonian is called the cluster Hamiltonian. 
The cluster Hamiltonian is given by the same form as 
Eq. 0 and now the indices i,j run only i,j = 
where Lc is the size of the finite cluster. When we in¬ 
vestigate the spontaneous symmetry breaking we include 
appropriate Weiss field terms into this cluster Hamilto¬ 
nian to resolve the degeneracy of the symmetry breaking 
directions. As an optional extension, Lf, sites are added 
at the boundary of the Lc-site cluster, which are called 
as “bath” sites and play a role to simulate the environ¬ 
ments. When the bath sites are added, we set up the 
cluster Hamiltonian on the (Lc + Lf,)-site cluster and in¬ 
clude new hybridization hopping terms so that electrons 
can move between the original Lc-site cluster and bath 
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sites. The electron correlations are neglected on these 
bath sites, therefore the on-site Coulomb interactions are 
absent on the bath sites in the cluster Hamiltonian. Then 
we approximate the electron correlations of the infinite 
system by the electron correlations among the Lc corre¬ 
lated sites. To be precise, we approximate the self-energy 
of the infinite system by the exact self-energy of the 
sites of the cluster Hamiltonian defined on the Lh) 
cluster. The additional elements in the cluster Hamil¬ 
tonian such as, parameters in the Weiss field terms and 
hybridization hopping between the correlated and bath 
sites, are determined based on self-consistency conditions 
or some variational principle. In principle the results of 
the above analysis converges to the exact results of the in¬ 
finite system in the large limit with Li,!Lf. —^ 0. Here, 

the latter condition is necessary because there are no fic¬ 
titious bath sites without electron correlations in the real 
system. In practice, the results depend on the cluster size 
and information on the infinite system is extracted using 
the cluster size dependence. In the present study, we need 
to compare the critical interaction strength Uaf of the 
magnetic transition and that of a non-magnetic metal-to- 
insulator transition Umi- In general, Uaf increases as the 
cluster size because the spatial fluctuations, which desta¬ 
bilize an ordered state, are taken into account more on a 
larger cluster. 17 mi also increases as the cluster size. This 
is because the electron wave functions can spread wider 
region on a larger cluster, which lowers the kinetic ener¬ 
gies and stabilizes a metallic state. Therefore we need to 
compute Uay and Umi at least on two clusters of different 
size and analyze their cluster size dependence. 

Next we briefly discuss the results of the previous in¬ 
vestigations on this matter and explain the controver¬ 
sial situations. Some years ago, Meng et al. (U found 
that Uaf = 4.3 and C/mi = 3.5 in model Q using 
the Quantum Monte Carlo simulations (QMC). Below 
U « 3.5, the system is a semi-metallic state, and be¬ 
yond U « 4.3 it is an antiferromagnetic state, and in the 
range 3.5 < t/ < 4.3 a spin liquid state is realized. Their 
results inspired a lot of theoretical interests. Sorella et 
al. [3] performed QMC for a system larger than Meng 
et al. They found that Uaf = 3.8 and the system is in 
a metallic state below Uaf- Thus their result ruled out 
the existence of a spin liquid state. The two subsequent 
QMC studies mis] support the result of Sorella et al. 

The quantum cluster methods are also applied to this 
model. Yu et oC[6] used VC A with the six-site hexago¬ 
nal cluster and reported that a spin liquid is realized for 
3.0 < C/ < 4.0. Wu et aC[7] used the cluster dynamical 
mean field theory (CDMFT) with clusters up to twenty- 
four sites and obtained e.g., Uaf = 3.7 and Umi = 3.2 
for the six-site hexagonal cluster, and Uaf = 3.8 and 
C^Mi = 3.2 for a twenty-four site cluster. Both Yu et 
al. and Wu et al. argued that their results are in good 
agreement with Meng et al. Seki and Ohta[S] studied 
this model using VC A and CDMFT and found that an 


insulating gap opens for infinitesimally small U on some 
clusters while Uaf remains finite. So their results also 
support the existence of a spin liquid state. However, 
this gap which opens at infinitesimally small U was fur¬ 
ther studied in Refs. [Ml] and the accumulated results 
of these studies indicate that this gap does not corre¬ 
spond to the gap due to the electron correlations, but 
it is an unphysical artifact arising because the trans¬ 
lational invariance is partially violated in the quantum 
cluster methods. Even though the precise condition is 
not yet clarified under which this unphysical gap opens 
at infinitesimally small U, so far this problem is observed 
only for some clusters in bipartite lattices. This unphysi¬ 
cal gap pushes electrons away from Fermi level and makes 
the formation of an ordered state difficult. For example, 
a nesting at Fermi level is lost. Therefore we can not 
study [/af,mi on a cluster with this problem. We need 
to analyze the cluster size dependence towards the limit 
Lc ^ oo with Lb/Lc —>■ 0 selecting the clusters which do 
not have this pathological problem. 

Hassan and SenechaljS] found that this unphysical gap 
is absent for the clusters {Lc,Lb) = (2,4) and (4,6), 
and analyzed [/af,im using CDMFT and cluster dynam¬ 
ical impurity approximation(CDIA) on these clusters. 
They found that Uaf — 3.3 and Umi — 5.5 for the 
{Lc,Lb) = (2,4) cluster, and Uaf — 4.0 and Umi — 6.3 
for the {Lc, Lb) = (4,6) cluster in CDMFT. In the case 
of CDIA, the metal-to-insulator transition becomes a 
first order transition and the coexisting range of the 
metallic and insulating phases are 5.7 ^ U ^ 6.0 for 
{Lc,Lb) = (2,4), and 6.5 < D < 8.0 for {Lc,Lb) = (4,6). 
Based on these results and observing that Umi is much 
larger than Uaf they ruled out the existence of a spin liq¬ 
uid phase. However, their analysis always uses the clus¬ 
ters satisfying Lb < Lc and largely violates the condition 
Lc > Lb. So it seems to be highly non-trivial if their 
results evaluate and approach towards the limit Lc ^ oo 
with LbjLc —>■ 0. In fact, quantitatively, their Umi on the 
{Lc, Lb) = (4,6) cluster is already much larger than the 
results of Meng et al. and Wu et al., and still increasing 
as the cluster size. 

In the case of VCA, Uaf — 1.5|8] and 17 mi — 3.4[TT] 
for the {Lc, Lb) = (8, 0) cluster, and 17 af — 2.7(5] and 
C7mi — 3.0[TT] for the {Lc, Lb) = (10, 0) cluster. 17mi de¬ 
creases rather largely as the cluster size increases, which 
is an anomalous behavior contradicting with the gen¬ 
eral argument of the cluster size dependence. Ebato et 
aL[T2| found that a gap opens at infinitesimally small U 
for {Lc,Lb) = (12,0) cluster in VCA. Therefore proper 
information on the cluster size dependence of Uaf,mi 
within VAC is still missing. We study the metal-to- 
insulator transition and magnetic transition by VCA us¬ 
ing {Lc,Lb) = (10,0) and (16,0) clusters in Fig. We 
found that Uaf = 2.7 and Umi = 3.0 for the 10-site clus¬ 
ter, and Uaf = 2.7 and Umi = 3.2 for the 16-site cluster. 
These results follow the general arguments of the cluster 
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size dependence of the C/af,mi and C/mi — Uaf increases 
slightly as the cluster size. This result rules out the ex¬ 
istence of a spin liquid state. 

The formulation .— In VGA we use the thermodynamic 
grand potential fit written in the form of a functional of 
the self-energy S as 

fitp] = Fm + Trln(-(Go 1 - E)-i) (2) 

for the Hamiltonian H of the infinite system and for the 
sum of the cluster Hamiltonians defined on each clus¬ 
ter tiling the infinite lattice, which we denote as H'. 
In Eq. (|^, Go is the non-interacting Green’s function, 
F[S] is the Legendre transform of the Luttinger-Ward 
functional [T5], and the index t denotes the explicit depen¬ 
dence of Ht on all the one-body operators in the Hamilto¬ 
nian. This functional is stationary (5r2t[S]/^S] = 0 at the 
true self-energy of the corresponding Hamiltonian and 
this stationary condition is equivalent to the Dyson’s 
equation. Bath sites are not used in VGA. All Hamil¬ 
tonians with the same interaction part share the same 
functional form of T'[E][T5]. Since the interaction part, 
which is the on-site Coulomb repulsion, is the same for 
H and H', F[E] is the same for a given E for FI and FI'. 
So eliminating /^[E], we obtain 

Dt[E] = D;,[E] + Trln(-(Go 1 - E)-i) 

-Trln(-(G'-i-E)-i), (3) 

where the prime denotes the quantities of FI'. This is 
a functional relation between Dt [E] and DJ., [E]. By ex¬ 
actly solving iL', we compute the Dt[E] for the exact 
self-energy of F['. To investigate symmetry breaking, we 
include the Weiss terms into F[' and the parameters in 
the Weiss terms are determined by solving the stationary 
condition. The exact self-energy of FI' satisfying the sta¬ 
tionary condition, denoted as E*, and r2t[E*] are the ap¬ 
proximate self-energy and grand-potential of FI in VGA. 
Physical quantities, such as expectation values of one- 
body operators, are evaluated using the Green’s function 
Go“^ — E*. In VGA, the short-range correlations within 
the cluster are exactly taken into account and the restric¬ 
tion of the space of the self-energies E into that of H' is 
the only approximation. 

In the present case the Weiss held is given by Haf = 
^AFl]iSign(i)(ni^ - rii^) with sign(z) = - 1 ( 1 ) for the 
hlled (unhlled) sites in Fig. We set the chemical po¬ 
tential of the system and cluster as /i = /i' = 17/2, 
which ensures the stationary half-hlling solutions be¬ 
cause the lattice is bipartite and particle-hole symme¬ 
try exists. The stationary condition is then reduced to 
dQ{hAF)/dhAF = 0 . In general, a non-magnetic solution 
with hAF = 0, and a magnetic solution with Haf 0 
are obtained and we compare their energies to determine 
which is the stable ground state. In our study it is suffice 
to compare the grand potential per site ^(/iaf), because 
/r is hxed to be t//2 at half-hlling. To further investigate 
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FIG. 2. (Color online) (a)-(b) The grand potential per 
site H(/iaf) as a function of /iaf with U = 2.7 (circles) and 
17 = 2.8 calculated on the 16-site-cluster. H(— /iaf) = H(/iaf) 
so only 0 < /iaf is shown, (c)-(f) The density of states 
around the Fermi level w = 0 computed on the 16-site cluster. 


the properties of the ground state, we also compute the 
density of state per site 

D[uj) = lim f ^ {--lmGaa{k,uj + iTj)X^) 

V->-0 / ZTT)^ rir ' TT 
^ ' a.a—l 

and the magnetic order parameter per site 

^ "c 

M = — EE sign(a)(ua^ - ria;), (5) 

a=l i 

where = 2 is the number of the sites in the unit cell in 
the sense of the sub-lattice formalism, and the k integra¬ 
tion is over the corresponding Brillouin zone. In Eq. 

77 —>■ 0 limit is evaluated using the standard extrapolation 
method by calculating D{uj) for rj = 0.01, 77 / 2 , and ry/4. 

Results. — Fig. (a)-(b) show the grand potential per 
site H(/iaf) as a function of Haf at (a) U = 2.7 and (b) 
U = 2.8 calculated on the 16-site cluster. A magnetic so¬ 
lution with /iaf 7 ^ 0 appears at U = 2.8 and the energy 
of this magnetic solution is lower than the non-magnetic 
solution with /iaf = 0 , so the magnetic state is realized 
at U = 2.8. For U < 2.7, a magnetic solution is not 
obtained and a non-magnetic state is realized. Fig. i(c)- 
(f) show the density of state D{uj) around the Fermi level 
a; = 0 as a function of oj calculated on the 16-site clus¬ 
ter imposing Haf = 0. For U < 3.2 the density of state 
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FIG. 3. (Color online) The phase diagram computed on 
the 16-site clnster. The filled circles are the density of state 
D{u} = 0.01) and filled triangles are the magnetic order 
parameter M. The unfilled marks correspond to the data 
computed on the 10-site cluster. 
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C^AF.Mi and [/mi — Uaf slightly increases as the cluster 
size. Therefore in this system, Umi and Uaf are rather 
close and still Uaf < Umi- The magnetic order parame¬ 
ter is not large but non-negligible in the region U ~ Umi- 
Thus a spin liquid state or a state very close to that is 
absent in this system. The analysis of the spectral weight 
function shows that the semi-metallic nature due to the 
Dirac cones are maintained up to the magnetic transition. 
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versity, ISSP, and Yukawa Institute. 


is very close to zero only at w = 0 around Fermi level 
and a semi-metallic nature is observed. For U > 3.3 a 
gap opens at Fermi level and the system is an insulator. 
By the similar analysis we found that the magnetic states 
are insulators. Fig.j^shows the phase diagram computed 
on the 16-site cluster. The filled circles are the density 
of state D{uj = 0.01) and filled triangles are the mag¬ 
netic order parameter M. D{uj = 0) vanishes both in a 
semi-metal and insulator so D{uj = 0.01) is shown to ob¬ 
serve an insulating gap. The unfilled marks correspond 
to the data computed on the 10 -site cluster. Uaf = 2.7 
and [/mi = 3.2 on the 16-site cluster and Uaf = 2.7 
and Umi = 3.0 on the 10 -site cluster. Our results of the 
10 -site cluster agree with the previous studies [ 8 l fTTl IT^ . 
Our results follow the general arguments on the cluster 
size dependence and in this sense anomalous behavior is 
not observed. The values of the order parameter sug¬ 
gests that the true minimum of the effective potential 
is slightly better simulated on the 16-site cluster. Both 
the magnetic and non-magnetic metal-to-insulator tran¬ 
sitions are of the second order since there are no coexist¬ 
ing regions. The magnetic order parameter is M ~ 1.4 
around U ~ Umi- 

Summary and discussion. — In conclusion, we have 
studied the possible emergence of a spin liquid state in 
the Hubbard model on the honeycomb lattice by VGA. 
We computed the critical interaction strength Uaf of a 
magnetic transition and that of the non-magnetic metal- 
to-insulator transition [/mi, and found that Uaf = 2.7 
and Umi = 3.2 on the 16-site cluster, and Uaf = 2.7 and 
Umi = 3.0 on the 10-site cluster. These results follow the 
general arguments of the cluster size dependence of the 
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